#include "NewtonInterp.h"

NewtonInterp::NewtonInterp(InterConditions p)
{
  diff_table=p.output_difftable();
  X_=p.get_X();
}

void NewtonInterp::Interpolation_Method1()
{
  vector<double> arr(1,diff_table[0][0]);
  Polynomial p(0,arr);
  vector<double> arr0(2,1);
  interPoly_=p;
  arr[0]=1;
  Polynomial p0(0,arr);
  for(int i=0;i<X_.size()-1;i++)
    {
      if(X_[i]==0) arr0[0]=X_[i];
      else  arr0[0]=-X_[i];
      Polynomial p1(1,arr0);
      p0=p0*p1;
      interPoly_=interPoly_+(diff_table[i+1][i+1])*p0;
    }
}

void NewtonInterp::Interpolation_Method2(vector<double> New_X,vector<double>New_f,vector<double> New_m,vector<vector<double> > New_df)
{
  int s=X_.size();
  int t=New_X.size()+New_m.size();
  diff_table.resize(s+t);
  for(int i=s;i<s+t;i++)
    diff_table[i].resize(i+1);
  X_.resize(s+t);
  int k=s;
  int index;
  for(int i=0;i<New_X.size();i++)
    for(int j=0;j<=New_m[i];j++)
      {
	X_[k]=New_X[i];
	diff_table[k][0]=New_f[i];
	k++;
      }
  for(int i=0;i<s+t-1;i++)
    for(int j=s;j<s+t;j++)
      {
        if(X_[j]!=X_[j-i-1])
  	  diff_table[j][i+1]=(diff_table[j][i]-diff_table[j-1][i])/(X_[j]-X_[j-i-1]);
       else
  	  {
	    for(int k=0;k<New_X.size();k++)
		if(X_[j]==New_X[k])
		  {break;index=k;}
 	    diff_table[j][i+1]=New_df[index][i]/factorial(i+1);
  	  }
      }

  vector<double> arr(1,1);
  vector<double> arr0(2,1);
  Polynomial p0(0,arr);
  for(int i=0;i<s-1;i++)
    {
      if(X_[i]==0) arr0[0]=X_[i];
      else  arr0[0]=-X_[i];
      Polynomial p1(1,arr0);
      p0=p0*p1;
    }
  for(int i=s-1;i<s+t-1;i++)
    {
      if(X_[i]==0) arr0[0]=X_[i];
      else  arr0[0]=-X_[i];
      Polynomial p1(1,arr0);
      p0=p0*p1;
      interPoly_=interPoly_+(diff_table[i+1][i+1])*p0;
    }
}





Polynomial NewtonInterp::get_interPoly_()const
{
  return interPoly_;
}

int NewtonInterp::factorial(int i)const
{
  int sum=1;
  for(int j=1;j<=i;j++)
    sum=sum*j;
  return sum;
}

double NewtonInterp::Neville_Aitken(InterConditions P, double x)
{
  int n=P.get_f().size();
  vector<double> _f=P.get_f();
  vector<double> _x=P.get_X();
  for(int i=0;i<n-1;i++)
    for(int j=n-1;j>=i+1;j--)
      _f[j]=(_f[j]*(x-_x[j-i-1])-_f[j-1]*(x-_x[j]))/(_x[j]-_x[j-i-1]);
  return _f[n-1]; 
}
